# -*- coding: utf-8 -*-
"""
Created on Tue Jun 17 13:41:57 2014

@author: Maxim
"""
from scipy.interpolate import interp1d
import numpy as np



def get_damping_in_roll_deriv(AR):
    """ NOTE: data only for taper=1.0
    Fig 15-18
    """
    _AR = np.array([6.02953, 10.0343, 13.1943, 16.0118])
    _Clr = np.array([-0.486278, -0.579635, -0.622251, -0.652031])
    curve = interp1d(_AR,_Clr, 'cubic')
    return curve(AR)


def get_CldeltaK(AR, b2start, b2end):
    _AR = np.array([6.0,10,16])
    _b2 = np.array([0.,0.2, 0.4, 0.6, 0.8, 1.0])
    node1 = np.array([5.94E-02, 4.96E-02, 3.69E-02])
    node2 = np.array([2.20E-01, 1.92E-01, 1.65E-01])
    node3 = np.array([4.72E-01, 4.26E-01, 3.69E-01])
    node4 = np.array([7.77E-01, 6.82E-01, 5.86E-01])
    node5 = np.array([9.81E-01, 8.87E-01, 7.71E-01])
    curve1 = interp1d(_AR,node1,'quadratic')
    curve2 = interp1d(_AR,node2,'quadratic')
    curve3 = interp1d(_AR,node3,'quadratic')
    curve4 = interp1d(_AR,node4,'quadratic')
    curve5 = interp1d(_AR,node5,'quadratic')
    pts = np.zeros(len(_b2))
    pts[1] = curve1(AR)
    pts[2] = curve2(AR)
    pts[3] = curve3(AR)
    pts[4] = curve4(AR)
    pts[5] = curve5(AR)
    curve6 = interp1d(_b2,pts,'cubic')
    return curve6(b2end) - curve6(b2start)


def get_k(ailRatio):
    _ailToChord = np.array([0.000000, 0.044692, 0.137309, 0.230276, 0.343041])
    _k = np.array([0.000000, 0.204369, 0.407511, 0.555129, 0.698643])
    curve = interp1d(_ailToChord,_k,'cubic')
    return curve(ailRatio)


def get_rate_of_roll(pb2V,Vt,span):
    """ figure 15-20 """
    if Vt<50: Vt = 50
    if Vt>500: Vt = 500
    if span<9.144: span=9.144
    if span>45.72: span=45.72
    
    _Vt = np.array([50, 100, 150, 200, 250, 300, 350, 400, 450, 500])
    _b = np.array([0.061935978, 0.03163129, 0.020536211, 0.015403197, 0.012208839, 0.009921189, 0.008428104, 0.007199392, 0.006288865, 0.005532865])
    _k = np.array([1.081141929, 0.564159023, 0.374882493, 0.28731583, 0.232821663, 0.193795486, 0.168324207, 0.147362998, 0.131829854, 0.118932859])
    curve1 = interp1d(_Vt,_b,'cubic')
    curve2 = interp1d(_Vt,_k,'cubic')
    b = curve1(Vt)
    k = curve2(Vt)
    x = (pb2V-b)/k
    
    _b2 = np.array([20.08157373, 17.0771128, 14.88449353, 13.2640195, 11.86201066, 10.69158747, 9.83437785, 8.331426206, 7.277668177, 6.427893896, 5.71934013, 5.206484495, 4.769669521, 4.354403597, 4.048759081, 3.734579837])
    _k2 = np.array([381.5805654, 325.4123404, 284.4214486, 254.1267795, 227.9163046, 206.0353101, 190.0098254, 161.9122305, 142.2122845, 126.325803, 113.0794309, 103.4916242, 95.32539328, 87.56202063, 81.84801391, 75.97445095])
    _span = np.array([9.144, 10.668, 12.192, 13.716, 15.24, 16.764, 18.288, 21.336, 24.384, 27.432, 30.48, 33.528, 36.576, 39.624, 42.672, 45.72])
    curve3 = interp1d(_span,_b2,'cubic')
    curve4 = interp1d(_span,_k2,'cubic')
    b2 = curve3(span)
    k2 = curve4(span)
    rateOfRoll = b2+k2*x
    return rateOfRoll

def rate_of_roll(ac,deflection=20.0, velocity=30):
    """ calculates rate of roll deg/sec"""
    AR = ac.wing.aspectRatio
    Clp = get_damping_in_roll_deriv(AR)
    b2s = 2.*ac.wing.aileron.bStart/ac.wing.span
    b2e = 2.*ac.wing.aileron.bEnd/ac.wing.span
    ailChordRatio = ac.wing.aileron.avgChordRatio
    CldeltaK = get_CldeltaK(AR,b2s,b2e)
    k = get_k(ailChordRatio)
    pb2V = -deflection*k/57.3*CldeltaK/Clp
    span = ac.wing.span
    return get_rate_of_roll(pb2V, velocity, span)


def run_test1():
    import aircraft
    ac = aircraft.load('V200')
    print rate_of_roll(ac,15,30), 'deg/sec'

if __name__=="__main__":
    run_test1()